#include<stdio.h>
#include<Algorithm/mesh_bool.h>
#include<Mesh_IO/Mesh_IO.h>

#include <time.h>
#include <Math/LB_Math.h>
#include <Matrix/LB_Matrix.h>
#include<Algorithm/subdivision_of_polygon.h>
#include<Mesh/libcell_iterator.h>
#include<Arcroll.h>

#include<libcell_macros.h>
static inline  int point_is_in_triangle(double *p1,double * p2,double* p3,double *p);

int is_point_in_area(double* p, double**ps,int len,int**subs);
double compute_function_value_on_mesh(Node* n,double *p);


static void my_sscanf(char * str,int* i1,int* i2,int *i3)
{
    int len=0;*i1=-1;*i2=-1;*i3=-1;
    char** words=split(str,'/',&len);
    int len1=strlen(str),num=0;
    for(int i=0;i<len1;i++)
    {
        if(str[i]=='/')
        {
            num++;
        }
    } 
    int * indexs[3]={i1,i2,i3};
    if(num==2&&len==2)
    {
        indexs[1]=i3; 
    }  
    for(int i=0;i<len;i++)
    {
        sscanf(words[i],"%d",indexs[i]);
        free(words[i]);
    }
    SAFE_FREE(words);
}
typedef struct Triangle{
	double ps[3][3];
}Triangle;

Node* get_some_triangles_of_intersections_from_file(double** ps,int len,int**subs,char* filename)
{
#ifdef _WIN32
    FILE *infile;
    fopen_s(&infile,filename,"r");
#else
    FILE *infile=fopen(filename,"r");
#endif
  	if(!infile)
    {
        printf("can't open this file %s\n",filename);
    }
    Node* re=NULL; char str[400]={NULL};
    fseek(infile,0,SEEK_SET);
   	Int_RB_Tree* tree_v=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
   	int_rb_tree_init(tree_v); 
   	int id=0,id_v=-1,id_vt=-1,id_vn=-1;
    int id_vs[3];
    double* p1=NULL,*p2=NULL,*p3=NULL;
    while(!feof(infile))
    {
 		fgets(str,400,infile);
        //printf("%s",str);
        if(feof(infile))
        {
            break;
        }
       	if(str[0]=='#')
        {
            continue;
        } 
    	else if(str[0]=='v'&&str[1]==' ')
        {
            double *point=(double*)malloc(sizeof(double)*3);
            sscanf(&str[2],"%lf %lf %lf",&(point[0]),&(point[1]),&(point[2]));
			if(is_point_in_area(point,ps,len,subs))
			{
            	tree_v->insert(tree_v,id,point);
			}
            id++;
        }
        else if(str[0]=='f'&&str[1]==' ')
        {
        	int fv_rows=0;
        	char**words=split(&str[2],' ',&fv_rows);
         	for(int i=0;i<fv_rows;i++)
            {
            	id_v=-1;id_vt=-1;id_vn=-1;
            	my_sscanf(words[i],&id_v,&id_vt,&id_vn);
                id_v--;id_vt--;id_vn--;
                id_vs[i]=id_v;
                free(words[i]);
            }
            // p1=(double*)tree_v->find(tree_v,id_vs[0]);
            // p2=(double*)tree_v->find(tree_v,id_vs[1]);
            // p3=(double*)tree_v->find(tree_v,id_vs[2]);
            // if(p1!=NULL||p2!=NULL||p3!=NULL)
            // {
            // 	Triangle *tri=(Triangle*)malloc(sizeof(Triangle));


            // }
            SAFE_FREE(words);
        }
    }
    int_rb_tree_free(tree_v);
    return re;
}

Mesh* get_cut_mesh_from_some_points(float** data,int len )
{
	if(len==0||data==NULL)
	{
		return NULL;
	}
	Mesh* re=(Mesh*)malloc(sizeof(Mesh));
	Mesh_init(re);
	template_v** vs=(template_v**)malloc(sizeof(template_v*)*len*2);
	double p[3]={0};
	for(int i=0;i<2*len;i++ )
	{
		p[0]=data[i%len][0];p[1]=data[i%len][1];p[2]=(i/len)*200-100;	
		vs[i]= re->create_vertexv(re,p,3);
	}

	template_v* vs1[3]={NULL};
	for(int i=0;i<len;i++)
	{
		vs1[0]=vs[i];vs1[1]=vs[(i+1)%len];vs1[2]=vs[i+len];
		re->create_cellv(re,vs1,3);
		vs1[0]=vs[(i+1)%len];vs1[1]=vs[(i+1)%len+len];vs1[2]=vs[i+len];
		re->create_cellv(re,vs1,3);
	}	
	free(vs);
	return re;
}



void test()
{
	Mesh mesh;
	Mesh_init(&mesh);
	_ReadObj_(&mesh,"texture_YAY_0108.obj");
	Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree);
	Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree1);



	int_rb_tree_free(tree1);
	int_rb_tree_free(tree);
	Mesh_free(&mesh);
}

static inline  int point_is_in_triangle(double *p1,double * p2,double* p3,double *p)
{
	double t1[3]={p[0]-p1[0],p[1]-p1[1],p[2]-p1[2]};
	double t2[3]={p[0]-p2[0],p[1]-p2[1],p[2]-p2[2]};
	double t3[3]={p[0]-p3[0],p[1]-p3[1],p[2]-p3[2]};
	double* n1=out_product(t1,t2);
	double* n2=out_product(t2,t3);
	double* n3=out_product(t3,t1);
	int re=0;
	if(inner_product(n1,n2,3)>0&&inner_product(n1,n3,3)>0 )
	{
		re= 1;
	}
	SAFE_FREE(n1);SAFE_FREE(n2);SAFE_FREE(n3);
	return re;
}

double* get_bounding_box_of_one_triangle(double* p1,double*  p2,double* p3)
{
	double* ps[3]={p1,p2,p3};
	double * re=(double*)malloc(sizeof(double)*4);
	re[0]=p1[0];re[1]=p1[1];re[2]=p1[0];re[3]=p1[1];
	for(int i=0;i<3;i++)
	{
		for(int j=0;j<2;j++)
		{
			re[j]= (ps[i][j]<re[j]?(ps[i][j]):re[j]);
			re[j+2]= (ps[i][j]>re[j+2]?(ps[i][j]):re[j+2]);
		}
	}
	return re;
}
double* get_bounding_box_of_one_polygon(double**ps,int len)
{
	double 	* re=(double*)malloc(sizeof(double)*4);
	re[0]=ps[1][0];re[1]=ps[1][1];re[2]=ps[1][0];re[3]=ps[1][1];
	for(int i=0;i<len;i++)
	{
		for(int j=0;j<2;j++)
		{
			re[j]= (ps[i][j]<re[j]?(ps[i][j]):re[j]);
			re[j+2]= (ps[i][j]>re[j+2]?(ps[i][j]):re[j+2]);
		}
	}
	return re;
}
double**  pre_get_pixel_data(double *boundings,double**ps,int len,int ** subs)
{
	int size=150;
	double** re=(double**)malloc(sizeof(double*)*size);
	for(int i=0;i<size;i++)
	{
		re[i]=(double*)malloc(sizeof(double)*size);
		memset(re[i],0,sizeof(double)*size);
	}	
	double stepx=(boundings[2]-boundings[0])/((double)size);
	double stepy=(boundings[3]-boundings[1])/((double)size);

	double area=stepx*stepy,sum=0,p[3]={0,0,0};	
	double* boundings1=NULL;
	int minx=0,miny=0,maxx=0, maxy=0;	
	for(int i=0;i<len-2;i++ )
	{
		boundings1=get_bounding_box_of_one_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]]);
		minx=(int)((boundings1[0]-boundings[0])/(stepx));
		miny=(int)((boundings1[1]-boundings[1])/(stepy));
		maxx=(int)((boundings1[2]-boundings[0])/(stepx))+1;
		maxy=(int)((boundings1[3]-boundings[1])/(stepy))+1;
		minx=((minx<0)?0:(minx));
		miny=((miny<0)?0:(miny));
		maxx=((maxx>150)?(150):maxx);
		maxy=((maxy>150)?(150):maxy);
		// if(minx>maxx||miny>maxy)
		// {
		// 	continue;
		// }
		for(int k=minx;k<maxx;k++)
		{
			p[0]=boundings[0]+stepx*(k+0.5);
			for(int j=miny;j<maxy;j++)
			{
				p[1]=boundings[1]+stepy*(j+0.5);
				if(re[k][j]!=-1&&point_is_in_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]],p))
				{
					re[k][j]=-1;
				}
			}
		}
		free(boundings1);
	}
	return re;
}


double** get_pixel_data(double*boundings,double**data,Node* n)
{
	int size=150;
	double** re=(double**)malloc(sizeof(double*)*size);
	for(int i=0;i<size;i++)
	{
		re[i]=(double*)malloc(sizeof(double)*size);
		memset(re[i],0,sizeof(double)*size);
	}
	double stepx=(boundings[2]-boundings[0])/((double)size);
	double stepy=(boundings[3]-boundings[1])/((double)size);

	double area=stepx*stepy,sum=0,p[3]={0,0,0},np[3]={0,0,0},xs[3]={0,0,0};	
	double* boundings1=NULL;

	int minx=0,miny=0,maxx=0, maxy=0;	

	LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
	lb_matrix_init_double(mat,3,3);	
	double*mat_data=(double*)(mat->data);
	mat_data[0*3+2]=0; mat_data[1*3+2]=0; mat_data[2*3+2]=1.0;
	for(Node*nit=n;nit!=NULL;nit=(Node*)(nit->Next) )
	{
		template_c* c=(template_c*)(nit->value);
		for(int j=0;j<3;j++)
		{
			mat_data[j*3+0]=c->vertices[1]->point[j]-c->vertices[0]->point[j];
			mat_data[j*3+1]=c->vertices[2]->point[j]-c->vertices[0]->point[j];
		}	
		LB_Matrix* inv_mat=mat->inverse(mat);
		if(inv_mat==NULL)
		{
			continue;
		}
		double *inv_mat_data=(double*)(inv_mat->data);
		boundings1=get_bounding_box_of_one_triangle(c->vertices[0]->point,c->vertices[1]->point,c->vertices[2]->point);
		minx=(int)((boundings1[0]-boundings[0])/(stepx));
		miny=(int)((boundings1[1]-boundings[1])/(stepy));
		maxx=(int)((boundings1[2]-boundings[0])/(stepx))+1;
		maxy=(int)((boundings1[3]-boundings[1])/(stepy))+1;
		minx=((minx<0)?0:(minx));
		miny=((miny<0)?0:(miny));
		maxx=((maxx>150)?(150):maxx);
		maxy=((maxy>150)?(150):maxy);
		np[2]=p[2]-c->vertices[0]->point[2];
		for(int i=minx;i<maxx;i++)
		{
			p[0]=boundings[0]+stepx*(i+0.5);
			np[0]=p[0]-c->vertices[0]->point[0];
			for(int j=miny;j<maxy;j++)
			{
				if(data[i][j]==0)
				{
					continue;
				}
				p[1]=boundings[1]+stepy*(j+0.5);
				np[1]=p[1]-c->vertices[0]->point[1];
				//此处应该是对单个三角形
				xs[0]=xs[1]=xs[2]=0;
				for(int k=0;k<3;k++)
				{
					xs[0]+=np[k]*inv_mat_data[0*3+k];
					xs[1]+=np[k]*inv_mat_data[1*3+k];
				}
				if(xs[0]>=0&&xs[1]>=0&&xs[0]+xs[1]<=1)
				{
					double temp=xs[0]*c->vertices[1]->point[2]+ xs[1]*c->vertices[2]->point[2]+(1-xs[0]-xs[1])*c->vertices[0]->point[2];
					
					if(data[i][j]==-1)
					{						
						re[i][j]=temp;
						data[i][j]=-2;
					}
					else if(temp>re[i][j])
					{
						re[i][j]=temp;
					}
				}
			}
		}
		lb_matrix_free(inv_mat);
		free(boundings1);
	}
	lb_matrix_free(mat);
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			data[i][j]==-2?(data[i][j]=-1):(data[i][j]=data[i][j]);
		}
	}
	return re;
}

int is_point_in_area(double* p, double**ps,int len,int**subs)
{
	for(int i=0;i<len-2;i++)
	{
		if(point_is_in_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]],p))
		{
			return 1;
		}	
	}
	return 0;
}	
double** get_pixel_datan(double*boundings,double**data,Node* n)
{
	int size=150;
	double** re=(double**)malloc(sizeof(double*)*size);
	for(int i=0;i<size;i++)
	{
		re[i]=(double*)malloc(sizeof(double)*size);
		memset(re[i],0,sizeof(double)*size);
	}
	double stepx=(boundings[2]-boundings[0])/((double)size);
	double stepy=(boundings[3]-boundings[1])/((double)size);

	double area=stepx*stepy,sum=0,p[3]={0,0,0},np[3]={0,0,0},xs[3]={0,0,0};	
	double* boundings1=NULL;
	double temp=0;
	int minx=0,miny=0,maxx=0, maxy=0;	

	LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
	lb_matrix_init_double(mat,3,3);	
	double*mat_data=(double*)(mat->data);
	mat_data[0*3+2]=0; mat_data[1*3+2]=0; mat_data[2*3+2]=1.0;
	for(Node*nit=n;nit!=NULL;nit=(Node*)(nit->Next) )
	{
		double* tri_c=(double*)(nit->value);
		mat_data[0*3+0]=tri_c[1*3+0]-tri_c[0*3+0];mat_data[1*3+0]=tri_c[1*3+1]-tri_c[0*3+1];mat_data[2*3+0]=tri_c[1*3+2]-tri_c[0*3+2];
		mat_data[0*3+1]=tri_c[2*3+0]-tri_c[0*3+0];mat_data[1*3+1]= tri_c[2*3+1]-tri_c[0*3+1];mat_data[2*3+1]= tri_c[2*3+2]-tri_c[0*3+2];

		LB_Matrix* inv_mat=mat->inverse(mat);
		if(inv_mat==NULL)
		{
			continue;
		}
		double *inv_mat_data=(double*)(inv_mat->data);
		boundings1=get_bounding_box_of_one_triangle(&tri_c[0*3],&tri_c[1*3],&tri_c[2*3]);
		minx=(int)((boundings1[0]-boundings[0])/(stepx));
		miny=(int)((boundings1[1]-boundings[1])/(stepy));
		maxx=(int)((boundings1[2]-boundings[0])/(stepx))+1;
		maxy=(int)((boundings1[3]-boundings[1])/(stepy))+1;
		minx=((minx<0)?0:(minx));
		miny=((miny<0)?0:(miny));
		maxx=((maxx>150)?(150):maxx);
		maxy=((maxy>150)?(150):maxy);
		np[2]=p[2]-tri_c[0*3+2];
		for(int i=minx;i<maxx;i++)
		{
			p[0]=boundings[0]+stepx*(i+0.5);
			np[0]=p[0]-tri_c[0*3+0];
			for(int j=miny;j<maxy;j++)
			{
				if(data[i][j]==0)
				{
					continue;
				}
				p[1]=boundings[1]+stepy*(j+0.5);
				np[1]=p[1]-tri_c[0*3+1];
				//此处应该是对单个三角形
				//xs[0]=xs[1]=xs[2]=0;
				xs[0]=np[0]*inv_mat_data[0*3+0]+np[1]*inv_mat_data[0*3+1]+np[2]*inv_mat_data[0*3+2]; 
				xs[1]=np[0]*inv_mat_data[1*3+0]+np[1]*inv_mat_data[1*3+1]+np[2]*inv_mat_data[1*3+2];
				if(xs[0]>=0&&xs[1]>=0&&xs[0]+xs[1]<=1)
				{
					temp=xs[0]*tri_c[1*3+2]+ xs[1]*tri_c[2*3+2]+(1-xs[0]-xs[1])*tri_c[0*3+2];
					
					if(data[i][j]==-1)
					{						
						re[i][j]=temp;
						data[i][j]=-2;
					}
					else if(temp>re[i][j])
					{
						re[i][j]=temp;
					}
				}
			}
		}
		lb_matrix_free(inv_mat);
		free(boundings1);
	}
	lb_matrix_free(mat);
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			data[i][j]==-2?(data[i][j]=-1):(data[i][j]=data[i][j]);
		}
	}
	return re;
}
Node* get_some_triangles_of_intersections_from_mesh(double* boundings,double* m,int len_m)
{
	Node* re=NULL;
	int num=len_m/9;
	for(int i=0;i<num;i++)
	{
		for(int j=0;j<3;j++)
		{
			if( boundings[0]<=m[i*9+j*3+0]&& m[i*9+j*3+0]<=boundings[2]&&  boundings[1]<=m[i*9+j*3+1]&&m[i*9+j*3+1]<=boundings[3] )
			{
				re=node_overlying(re,&m[i*9]);
				break;
			}	
		}	
	}	
	return re;
}

Node* get_some_cells_of_intersections_from_mesh(double** ps,int len,int**subs,Mesh* m)
{
	Node* re=NULL;
	Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
	int_rb_tree_init(tree);
	double p[3]={0,0,0};
	for(auto vit=m->v_begin(m);vit!=m->v_end(m);vit++)
	{
		//printf("once ");
		p[0]=quote(vit)->point[0];p[1]=quote(vit)->point[1];
		if(is_point_in_area(p,ps,len,subs))
		{
			for(auto vcit=m->vc_begin(m,*vit);vcit!=m->vc_end(m,*vit);vcit++)
			{
				tree->insert(tree,quote(vcit)->id,quote(vcit));
			}
		}	
	}
	for(auto it=tree->begin(tree);it.it!=NULL;it++)
	{
		re=node_overlying(re,it.second);
	}
	int_rb_tree_free(tree);
	return re;
}




double compute_function_value_on_mesh(Node* n,double *p)
{
	double re=0;
	double np[3]={0};
	LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
	lb_matrix_init_double(mat,3,3);
	double * data=(double*)(mat->data),*data1=NULL;
	data[0*3+2]=0; data[1*3+2]=0; data[2*3+2]=1.0;
	double xs[3]={0} ;			
	for(Node *nit=n;nit!=NULL;nit=(Node*)(nit->Next))
	{
		template_c* c=(template_c*)(nit->value);
		for(int j=0;j<3;j++)
		{
			data[j*3+0]=c->vertices[1]->point[j]-c->vertices[0]->point[j];
			data[j*3+1]=c->vertices[2]->point[j]-c->vertices[0]->point[j];
			np[j]=p[j]-c->vertices[0]->point[j];
		}
		LB_Matrix* inv=mat->inverse(mat);
		data1=(double*)(inv->data);
		for(int i=0;i<3;i++)
		{	
			xs[i]=0;
			for(int j=0;j<3;j++)
			{	
				xs[i]+=np[j]*data1[i*3+j];
			}
		}	
		lb_matrix_free(inv);
		if(xs[0]>=0&&xs[1]>=0&&xs[0]+xs[1]<=1)
		{
			re= xs[0]*c->vertices[1]->point[2]+ xs[1]*c->vertices[2]->point[2]+(1-xs[0]-xs[1])*c->vertices[0]->point[2];
			break;
		}
	}
	lb_matrix_free(mat);
	return re;
}

double compute_integral_on_one_triangle(double* p1,double*p2,double*p3,Node* n)
{
	double * boundings=get_bounding_box_of_one_triangle(p1,p2,p3);
	printf("boundings:%lf %lf %lf %lf \n",boundings[0],boundings[1],boundings[2],boundings[3]);

	int fens=150;
	double stepx=(boundings[2]-boundings[0])/((double)fens);
	double stepy=(boundings[3]-boundings[1])/((double)fens);

	double area=stepx*stepy,sum=0,p[3]={0,0,0};
	for(int i=0;i<fens;i++)
	{
		p[0]=boundings[0]+stepx*(i+0.5);
		for(int j=0;j<fens;j++)
		{
			p[1]=boundings[1]+stepy*(j+0.5);

			sum+=compute_function_value_on_mesh(n,p);
		}
	}
	free(boundings);
	return sum*area;
}
void test2(double**ps,int len,char* filename1,char*filename2)
{
	int**subs=(int**)malloc(sizeof(int*)*(len-2));
	for(int i=0;i<len-2;i++)
	{
		subs[i]=(int*)malloc(sizeof(int)*3);
		memset(subs[i],0,sizeof(int)*3);
	}
	optimization_subdivision_of_polygon(ps,len,subs);
	printf("begin\n");

	Node* n1=get_some_triangles_of_intersections_from_file(ps,len,subs,filename1);
	Node* n2=get_some_triangles_of_intersections_from_file(ps,len,subs,filename2);


	free_node(n1);free_node(n2);
	for(int i=0;i<len-2;i++)
	{
		free(subs[i]);
	}
	free(subs);
}
static  void test_show_mesh_triangles(Viewer_World* vw,Node* node)
{
 	Node* n=vw->create_something(vw,"Faces");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
     
    vf->Data_rows=node_size(node)*3;
    vf->Data_index_rows=node_size(node);
   	vf->Data= (float*)malloc(sizeof(float)*vf->Data_rows*3);
   	vf->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*vf->Data_index_rows*4);
   	int i=0;
   	for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
   	{
   		double* tri_c=(double*)(nit->value);
   		for(int j=0;j<3;j++)
   		{
   			vf->Data[i*9+j*3+0]=tri_c[j*3+0];
			vf->Data[i*9+j*3+1]=tri_c[j*3+1];
   			vf->Data[i*9+j*3+2]=tri_c[j*3+2];
   		}	
   		vf->Data_index[i*4+0]=3;vf->Data_index[i*4+1]=i*3+0;vf->Data_index[i*4+2]=i*3+1;vf->Data_index[i*4+3]=i*3+2;
   		i++;
   	}
    vf->color_rows=vf->Data_index_rows;
    vf->set_color(vf,0.9,0.0,0.0,1.0);

    vf->normal_rows=vf->Data_index_rows;
    
    free_node(n);
}

static void  test_show_mesh_cells(Viewer_World* vw,Node* node)
{
    Node* n=vw->create_something(vw,"Faces");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
     
    vf->Data_rows=node_size(node)*3;
    vf->Data_index_rows=node_size(node);
   	vf->Data= (float*)malloc(sizeof(float)*vf->Data_rows*3);
   	vf->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*vf->Data_index_rows*4);
   	int i=0;
   	for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
   	{
   		template_c* c=(template_c*)(nit->value);
   		for(int j=0;j<3;j++)
   		{
   			vf->Data[i*9+j*3+0]=c->vertices[j]->point[0];
			vf->Data[i*9+j*3+1]=c->vertices[j]->point[1];
   			vf->Data[i*9+j*3+2]=c->vertices[j]->point[2];
   		}	
   		vf->Data_index[i*4+0]=3;vf->Data_index[i*4+1]=i*3+0;vf->Data_index[i*4+2]=i*3+1;vf->Data_index[i*4+3]=i*3+2;
   		i++;
   	}
    vf->color_rows=vf->Data_index_rows;
    vf->set_color(vf,0.9,0.0,0.0,1.0);

    vf->normal_rows=vf->Data_index_rows;
    
    free_node(n); 
}
static void test_show_intergral_points(Viewer_World*vw ,double**data,double*boundings)
{
	int size=150;
	double stepx=(boundings[2]-boundings[0])/((double)size);
	double stepy=(boundings[3]-boundings[1])/((double)size);
	
	Node* n=vw->create_something(vw,"Points");
	Viewer_Something* vs= (Viewer_Something*)(n->value);
	Viewer_Points* vp=(Viewer_Points*)(vs->evolution);

	vp->set_data_rows(vp,size*size);
	int num=0;	
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			if(data[i][j]<0)
			{
				vp->Data[num*3+0]=boundings[0]+stepx*(i+0.5);vp->Data[num*3+1]=boundings[1]+stepy*(j+0.5);
				//vp->Data[num*3+2]=data[i][j];	
				num++;
			}
		}
	}
	vp->Data_rows=num;	
	vp->set_color(vp,0.1,0.7,0.9,1.0);
	free_node(n);	
}
static void test_show_intergral_points1(Viewer_World*vw ,double**data,double*boundings)
{
	int size=150;
	double stepx=(boundings[2]-boundings[0])/((double)size);
	double stepy=(boundings[3]-boundings[1])/((double)size);
	
	Node* n=vw->create_something(vw,"Points");
	Viewer_Something* vs= (Viewer_Something*)(n->value);
	Viewer_Points* vp=(Viewer_Points*)(vs->evolution);
	vp->set_data_rows(vp,size*size);
	int num=0;
	for(int i=0;i<size;i++)
	{
		for(int j=0;j<size;j++)
		{
			vp->Data[num*3+0]=boundings[0]+stepx*(i+0.5);vp->Data[num*3+1]=boundings[1]+stepy*(j+0.5);
			vp->Data[num*3+2]=data[i][j];	
		
			num++;
		}
	}
	vp->Data_rows=num;	
	vp->set_color(vp,0.1,0.7,0.9,1.0);
	free_node(n);	
}
double*  my_test2( double* m1,int len_m1,double**ps ,int len )
{
	printf("begin\n");
	int**subs=(int**)malloc(sizeof(int*)*(len-2));
	for(int i=0;i<len-2;i++)
	{
		subs[i]=(int*)malloc(sizeof(int)*3);
		memset(subs[i],0,sizeof(int)*3);
	}
	optimization_subdivision_of_polygon(ps,len,subs);


	double* boundings =get_bounding_box_of_one_polygon(ps,len);
	//Node* n1=get_some_cells_of_intersections_from_mesh(ps,len,subs,m1);
    Node*n1=get_some_triangles_of_intersections_from_mesh(boundings,m1,len_m1);
	double**pixel_data1= pre_get_pixel_data(boundings,ps,len,subs);
	

	double**int_pixel_data1= get_pixel_datan(boundings,pixel_data1,n1);

	Viewer_World_Manager vwm;
	viewer_world_manager_init(&vwm);

	Viewer_World* vw=vwm.create_world(&vwm,NULL);
	add_default_somethings(vw);

	Node* n=vw->create_something(vw,"Edges");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
    ve->Data_rows=len;
    ve->Data_index_rows=len;
    ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
   	ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows); 
   	for( int i=0;i<len;i++)
   	{
   		ve->Data[i*3+0]=ps[i][0];ve->Data[i*3+1]=ps[i][1];	
   		ve->Data[i*3+2]=0; //ve->Data[i*3+2]=-52.0;
   		ve->Data_index[i*2+0]=i;ve->Data_index[i*2+1]=(i+1)%len;
   	}

    //get_lines_data_from_2dim_cell(m,&(ve->Data),&(ve->Data_index));
    ve->color_rows=ve->Data_index_rows;
    ve->set_color(ve,1.0,0.0,0.6,1.0); 
    free_node(n);

	//test_show_intergral_points(vw,pixel_data1 ,boundings );
	test_show_intergral_points1(vw,int_pixel_data1,boundings);
	test_show_mesh_triangles(vw,n1);


	Viewer_Opengl_Interpreter voi;
	viewer_opengl_interpreter_initn(&voi,&vwm);
	voi.interpreter(&voi);
	viewer_world_manager_free(&vwm);


	int size=150;
	double * re=(double*)malloc(sizeof(double)*(3*size*size+2));
	memset(re,0,sizeof(double)*(3*size*size+2));
	double stepx=(boundings[2]-boundings[0])/((double)size);
	double stepy=(boundings[3]-boundings[1])/((double)size);
	double area=stepx*stepy,temp=0,sum=0;int temp_num=0;
	for(int i=0;i<size;i++)
	{
		temp=(i+0.5)*stepx+boundings[0];
		for(int j=0;j<size;j++)
		{
			if(pixel_data1[i][j]==-1)
			{
				re[temp_num*3+0 +2]=temp;
				re[temp_num*3+1 +2]=(j+0.5)*stepy+boundings[1];
				re[temp_num*3+2+2]=int_pixel_data1[i][j];
				sum+=int_pixel_data1[i][j];
				temp_num++;
			}	
		}
		SAFE_FREE(pixel_data1[i]);
		SAFE_FREE(int_pixel_data1[i]);
	}
	re[0]=area*sum;
	re[1]=temp_num;
	//printf("cut area: %lf %lf\n",temp_num*area, compute_area_of_triangle(ps[0],ps[1],ps[2])+ compute_area_of_triangle(ps[2],ps[3],ps[0]));


	SAFE_FREE(int_pixel_data1);
	SAFE_FREE(pixel_data1);
	free_node(n1);
	for(int i=0;i<len-2;i++)
	{
		free(subs[i]);
	}
	free(subs);
	return re;
}

double*  my_test1( double* m1,int len_m1,double* m2,int len_m2,double**ps ,int len )
{
	double * re=NULL;
	printf("begin\n");
	int**subs=(int**)malloc(sizeof(int*)*(len-2));
	for(int i=0;i<len-2;i++)
	{
		subs[i]=(int*)malloc(sizeof(int)*3);
		memset(subs[i],0,sizeof(int)*3);
	}
	optimization_subdivision_of_polygon(ps,len,subs);


	double* boundings =get_bounding_box_of_one_polygon(ps,len);
	//Node* n1=get_some_cells_of_intersections_from_mesh(ps,len,subs,m1);
    Node*n1=get_some_triangles_of_intersections_from_mesh(boundings,m1,len_m1);
	Node*n2=get_some_triangles_of_intersections_from_mesh(boundings,m2,len_m2);
	double**pixel_data1= pre_get_pixel_data(boundings,ps,len,subs);

	double**int_pixel_data1= get_pixel_datan(boundings,pixel_data1,n1);
	double**int_pixel_data2= get_pixel_datan(boundings,pixel_data1,n2);






		

	SAFE_FREE(int_pixel_data1);
	SAFE_FREE(int_pixel_data2);
	SAFE_FREE(pixel_data1);
	free_node(n1);
	free_node(n2);
	for(int i=0;i<len-2;i++)
	{
		free(subs[i]);
	}
	free(subs);
	return re;
}


void test1(double**ps,int len,Mesh*m1,Mesh* m2)
{
	printf("begin\n");

	int**subs=(int**)malloc(sizeof(int*)*(len-2));
	for(int i=0;i<len-2;i++)
	{
		subs[i]=(int*)malloc(sizeof(int)*3);
		memset(subs[i],0,sizeof(int)*3);
	}
	optimization_subdivision_of_polygon(ps,len,subs);


	double* boundings =get_bounding_box_of_one_polygon(ps,len);
	//求n1 n2，可以根据boundings加速粗略判断
	Node* n1=get_some_cells_of_intersections_from_mesh(ps,len,subs,m1);
	Node* n2=get_some_cells_of_intersections_from_mesh(ps,len,subs,m2);
	//求n1 n2，可以根据boundings加速粗略判断
	printf("boundings:%lf %lf %lf %lf\n",boundings[0],boundings[1],boundings[2],boundings[3]);

	double**pixel_data1= pre_get_pixel_data(boundings,ps,len,subs);

	double**int_pixel_data1= get_pixel_data(boundings,pixel_data1,n1);
	double**int_pixel_data2= get_pixel_data(boundings,pixel_data1,n2);


	printf("n1 size:%d  n2 size:%d\n",node_size(n1),node_size(n2));

	Viewer_World_Manager vwm;
	viewer_world_manager_init(&vwm);

	Viewer_World* vw=vwm.create_world(&vwm,NULL);
	add_default_somethings(vw);
	test_show_mesh_cells(vw,n1);
	test_show_mesh_cells(vw,n2);
	//test_show_intergral_points(vw,pixel_data1 ,boundings );
	test_show_intergral_points1(vw,int_pixel_data2,boundings);

	Node* n=vw->create_something(vw,"Edges");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
    ve->Data_rows=len;
    ve->Data_index_rows=len;
    ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
   	ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows); 
   	for( int i=0;i<len;i++)
   	{
   		ve->Data[i*3+0]=ps[i][0];ve->Data[i*3+1]=ps[i][1];	
   		ve->Data[i*3+2]=0; //ve->Data[i*3+2]=-52.0;
   		ve->Data_index[i*2+0]=i;ve->Data_index[i*2+1]=(i+1)%len;
   	}

    //get_lines_data_from_2dim_cell(m,&(ve->Data),&(ve->Data_index));
    ve->color_rows=ve->Data_index_rows;
    ve->set_color(ve,1.0,0.0,0.6,1.0); 
    free_node(n);





	Viewer_Opengl_Interpreter voi;

	viewer_opengl_interpreter_initn(&voi,&vwm);
	voi.interpreter(&voi);

	viewer_world_manager_free(&vwm);

	double integral1=0,integral2=0;

	for(int i=0;i<len-2;i++)
	{
		//integral1+=compute_integral_on_one_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]],n1);
	// 	integral2+=compute_integral_on_one_triangle(ps[subs[i][0]],ps[subs[i][1]],ps[subs[i][2]],n2);
	}
	free_node(n1);
	free_node(n2);
	printf("%d %d %d\n",subs[0][0],subs[0][1],subs[0][2]);
	printf("%d %d %d\n",subs[1][0],subs[1][1],subs[1][2]);

	for(int i=0;i<len-2;i++)
	{
		free(subs[i]);
	}
	free(subs);
}


int main()
{
	//test();
 	
	double**ps=(double**)malloc(sizeof(double*)*4);
	for(int i=0;i<4;i++)
	{
		ps[i]=(double*)malloc(sizeof(double)*3);
		memset(ps[i],0,sizeof(double)*3);
	}
	ps[0][0]=240;ps[0][1]=40;ps[0][2]=0;
	ps[1][0]=300;ps[1][1]=40;ps[1][2]=0;
	ps[2][0]=290;ps[2][1]=60;ps[2][2]=0;
	ps[3][0]=260;ps[3][1]=70;ps[3][2]=0;

	Mesh mesh1,mesh2;
	Mesh_init(&mesh1);Mesh_init(&mesh2);
	_ReadObj_(&mesh1,"texture_YAY_0108.obj");
	_ReadObj_(&mesh2,"texture_YAY.obj");

	float a=1.99;
	int b=(int)a;
	printf("b:%d\n",b);
	double * array_m1=(double*)malloc(sizeof(double)*9*mesh1.num_c(&mesh1));
	int num_i=0;
	for(auto cit=mesh1.c_begin(&mesh1);cit!=mesh1.c_end(&mesh1);cit++)
	{
		for(int j=0;j<3;j++)
		{
			for(int k=0;k<3;k++)
			{
				array_m1[num_i*9+j*3+k]=quote(cit)->vertices[j]->point[k];
			}
		}
		num_i++;
	}	
	double * array_m2=(double*)malloc(sizeof(double)*9*mesh2.num_c(&mesh2));
	num_i=0;
	for(auto cit=mesh2.c_begin(&mesh2);cit!=mesh2.c_end(&mesh2);cit++)
	{
		for(int j=0;j<3;j++)
		{
			for(int k=0;k<3;k++)
			{
				array_m2[num_i*9+j*3+k]=quote(cit)->vertices[j]->point[k];
			}
		}
		num_i++;
	}	

	//my_test1( array_m1,9*mesh1.num_c(&mesh1),array_m2,9*mesh2.num_c(&mesh2),ps ,4 );
	double * re=my_test2(array_m1,9*mesh1.num_c(&mesh1),ps,4);


	


 	//test1(ps,4,&mesh1,&mesh2);

	//test2(ps,4,"texture_YAY.obj","texture_YAY_0108.obj");

	Mesh_free(&mesh1);Mesh_free(&mesh2);
	printf("end\n");
	return 0;
}
